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Abstract 

A  proton  exchange  membrane  fuel  cell  (PEMFC)  stack  model  which  incorporates  flow  distribution  effects  and  a  reduced  dimensional  unit 
cell  model  is  presented.  The  laws  of  Mass  and  Momentum  Conservation  are  applied  throughout  the  stack.  Along  the  headers,  flow  splitting  and 
recombination  are  applied  at  each  tee  junction,  while  along  the  unit  cell  channels,  reactant  consumption  and  byproduct  production  are  accounted 
for.  Cell  performance  is  coupled  to  the  channel  conditions.  Using  this  stack  model,  the  sensitivity  of  stack  performance  to  operating  conditions 
(inlet  velocity  and  pressure)  and  design  parameters  (manifold,  flow  configuration  and  friction  factor)  is  investigated.  In  particular,  performance 
under  uniform,  single  anomaly,  and  random  parameter  distributions  is  investigated. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Proton  exchange  membrane  fuel  cells  continue  to  attract  sig¬ 
nificant  attention  because  of  their  multiple  advantages,  including 
high  power  density,  low  pollution,  and  operation  at  relatively  low 
temperature  and  pressure.  A  significant  development  area  is  rep¬ 
resented  by  modeling  due  to  its  potential  to  reduce  development 
costs  and  accelerate  design  cycle  times.  This  is  especially  true 
for  commercial  sized  stacks,  where  a  significant  portion  of  de¬ 
velopment  costs  is  due  to  manufacturing,  assembly,  and  testing. 
However,  cell  interactions  within  stacks  have  not  yet  drawn  sig¬ 
nificant  attention,  especially  the  thermal  and  electrical  aspects 
[1,5,9,13],  yet  an  understanding  of  these  interactions  could  be 
key  towards  optimizing  stack  design  and  understanding  the  fun¬ 
damental  relationship  between  part  tolerances  and  performance. 
In  particular,  a  prime  design  objective  is  to  achieve  uniform  op¬ 
erating  conditions  across  all  cells,  as  stack  failure  is  dictated 
in  many  cases  by  the  weakest  cell  (series  reliability).  Evidence 
for  the  existence  of  non-uniform  cell  performance  distribution 
is  available  [12,14]. 

The  cell  interaction  which  has  received  the  most  attention  is 
flow  distribution  [3,4,7,8,10,1 1,17].  While  much  work  has  been 
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completed  on  the  effects  of  friction  loss  on  flow  distribution 
(including  wall  friction,  geometrical  friction  due  to  fittings  or 
valves,  etc.),  the  effects  of  reactant  consumption  in  the  unit  cell 
flow  field  channels  has  largely  been  ignored  [3,8, 1 1]  or  idealized 
[4,7].  For  instance,  one  idealization  is  the  assumption  that  the 
reactants  are  consumed  at  a  location  corresponding  to  half  the 
flow  field  length  [4,7].  More  realistic  unit  cell  models,  including 
[22] ,  have  not  yet  been  implemented  into  stack  environments,  in¬ 
cluding  features  such  as  consumption/generation  and  their  effect 
on  heat/mass  transfer.  Additionally,  the  effect  of  part  tolerances 
has  hardly  been  taken  into  account  [8]  but  was  recognized  as 
an  important  aspect  [5]  which  justifies  the  need  for  fast,  lower 
order  models  (1+1D  rather  than  3D)  to  balance  computational 
requirements  and  cell  monitoring  to  diagnose  issues.  Despite  the 
recognition  of  the  importance  of  manufacturing  tolerances,  pa¬ 
rameter  variations  were  not  taken  into  account  and  the  reactant 
flow  distribution  was  prescribed  rather  than  computed.  However, 
the  results  did  show  that  different  cell  coupling  effects  could  in¬ 
teract.  Finally,  the  flow  distribution  does  not  seem  to  have  been 
related  to  cell  performance,  at  least  for  the  cases  that  do  not  have 
a  specified  flow  distribution,  yet  this  is  a  requirement  for  under¬ 
standing  the  relationship  between  manufacturing  specifications 
and  product  performance. 

In  the  present  paper,  a  previous  unit  cell  performance  model 
including  consumption  of  reactants  and  product  generation  [2] 
is  coupled  to  a  stack-level  flow  distribution  model.  Particular 
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Nomenclature 

Symbol 

A 

cross-sectional  area  of  channel  or  header  (m2) 

a 

fixed  charge  (sulfonic  acid  group)  concentration 

Q 

(molm  ) 

C 

gas  concentration  (molm-3) 

Gref 

reference  oxygen  concentration  (molm-3) 

6‘w 

free  membrane  water  concentration  (molm-3) 

Ca 

membrane  water  content  at  anode  side  (mol  m-3) 

Cc 

membrane  water  content  at  cathode  side 

Q 

(molm  ) 

C+ 

membrane  hydronium  concentration  (molm-3) 

£>h 

hydraulic  diameter  (m) 

D+ 

hydronium  diffusion  coefficient  (m2  s-1) 

dw 

water  diffusion  coefficient  (m2  s-1) 

Eo 

open  circuit  voltage  (V) 

T 

Faraday’s  constant  (C  mol-1) 

f 

friction  factor  (— ) 

I 

current  (A) 

i 

current  density  (Am-2) 

io 

cathode  exchange  current  density  (Am-2) 

j 

index:  o  =  O2,  n  =  N2,  h  =  H2,  w  =  H2O 

7W 

water  crossover  flux  (molm-2  s-1) 

Gch 

length  of  channel  (m) 

Lh 

length  of  header  (m) 

L  m 

membrane  thickness  (m) 

G  w 

width  of  MEA  per  channel  (m) 

M 

molar  mass  (kg  mol-1) 

m 

mass  (kg) 

N 

number  of  unit  cells  (— ) 

number  of  channels  per  unit  cell  (— ) 

P 

pressure  (Pa) 

P 

cross-sectional  perimeter  of  channel  or  header 
(m) 

Q 

molar  flow  rate  (mol  s-1) 

R 

membrane  ohmic  resistance  (£2  •  m2) 

n 

Ideal  Gas  Constant  (J  mol-1  K-1) 

Re 

Reynold’s  number  (— ) 

s 

stoich  (— ) 

T 

temperature  (K) 

u 

gas  velocity  (/xs-1) 

V 

voltage  (V) 

y 

position  along  the  channel  (m) 

z 

position  along  the  header  (m) 

Symbol 

a 

water  crossover  from  anode  to  cathode  (— ) 

ac 

cathode  transfer  coefficient  (— ) 

S 

oxygen  mass  transfer  coefficient  (mol  A-1  m-1) 

P 

mass  density  (kgm-3) 

attention  is  given  to  model  flexibility  with  the  ability  to  vary 
all  specified  cell  parameters  independently  for  each  individual 
cell  and  potentially  reproduce  observed  manufacturing  varia¬ 
tions.  Computed  examples  are  discussed  focusing  on  the  oxidant 
stream  and  cell/stack  design  characteristics  related  to  flow.  Some 
validation  data  are  also  discussed.  Additional  effects  including 
two-phase  flow  in  the  channel  and  gas  diffusion  layer  (GDL), 
and  electrical  and  thermal  coupling  effects  are  not  considered 
in  this  present  model,  but  work  is  planned  to  incorporate  these 
effects  into  the  next  version  of  this  stack  model. 

2.  Mathematical  model 

2.7.  Unit  cell  model 

The  essential  features  of  a  previous  unit  cell  model  are  sum¬ 
marized  [2].  Oxygen  is  consumed  along  the  cathode  channel 
and  hydrogen  is  consumed  along  the  anode  channel  to  generate 
current.  The  associated  electrochemical  reactions  are 

*  4H+  +4e~ 


2H2  - 
02  +  4H+  +  4e 


2H20. 


(1) 

(2) 


respectively.  The  associated  mass  balance  equations  along  each 
channel  are 


d  Q 


o 


i(y)L 


W 


d y 

d<2h 
d  y 


=  =F 


AT 

i(y)Lw 


2  T 


(3) 


(4) 


where  the  minus  and  plus  signs  in  Eq.  (4)  are  taken  for  coflow 
and  counterflow  operation  respectively.  Here,  Lw  is  the  mem¬ 
brane/electrode  assembly  (MEA)  width  associated  with  each 
channel.  In  addition,  the  mass  balance  equations  for  water  in  the 
cathode  and  anode  channels  are 


d  Q 


w 


d y 

d<2w 

dy 


=  +(l+«O0) 


i(y)L' 

2T 


=  T  a(y) 


i(y)L. 
2  T 


(5) 


(6) 


respectively,  where  a  =  2TJW  /  i  is  a  measure  of  the  total  water 
crossover  flux  from  anode  to  cathode.  Jw  is  given  by  [2] 

(cl)1  -  (cl)1  '• 


J w  —  cidwj  e 


w 


-2436/  T 


2  Lm 

The  voltage  balance  is  [2] 

1ZT  _  iCVQ f 


+  T 


(7) 


V  =  E0-iR- 


ln 


(8) 


acT  io(C0  -  Si) 
where  the  membrane  protonic  resistance  is  given  by 

7 ZT  fLm  dy 

R=  — T  /  - - - .  (9) 

In  Eq.  (8),  the  cathode  overpotential  has  been  approximated  by 
its  Tafel  form  valid  for  high  currents.  The  term  C0  —  Si  approxi¬ 
mates  the  oxygen  concentration  at  the  catalyst  sites,  and  is  lower 
than  the  channel  average  because  oxygen  must  diffuse  through 
the  GDL  and  catalyst  layer  to  the  active  sites.  The  concentra¬ 
tion  reduction  should  be  proportional  to  the  flux,  which  in  turn 
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is  proportional  to  the  local  current,  leading  to  the  given  form. 
Eq.  (8)  holds  for  i  <  /max  =  C0/8,  where  /max  is  the  maximum 
current  allowed  locally  by  mass  transport. 

At  saturated  conditions,  the  channel  oxygen  concentration 
C0  is  related  to  the  oxygen  channel  molar  flow  rate  Go  by 


(P~P^t(T))Q0 

(Go  +  Qn)nr  ’ 


(10) 


where  it  is  assumed  that  the  Ideal  Gas  Law  applies,  that  oversat¬ 
urated  vapour  immediately  condenses  (no  liquid  water  effects 
are  accounted  for),  and  that  all  gaseous  species  are  transported 
by  channel  convection  with  a  common  average  velocity.  At  un¬ 
dersaturated  conditions, 


QqP 

(Go  H"  Gn  H"  Qw)KT 


(11) 


where  Gw  is  the  vapor  water  flux.  The  anode  gas  concentrations 
have  similar  expressions.  Note  that  Eqs.  (8)  and  (10)  couple 
the  cell  performance  to  the  channel  conditions.  The  saturation 
pressure  depends  on  temperature,  and  is  given  by  the  empirical 
formula  [15] 


Psat(T)  =  101325 

xlQ-2.i8+o.029(r-rc)-9.i8xio-5(r-rc)2+i.44xio-7(r-rc)3  ^2) 


where  Tc  is  the  ice  point  of  water. 

A  linear  temperature  profile  is  imposed  along  the  channel,  as 
well  as  Momentum  Conservation,  see  Eq.  (18).  The  current  den¬ 
sity  profile  must  satisfy  the  requirement  that  its  average  equals 
the  target  current  density 

t  rLch 

—  /  /(V)dv  =  /'avg-  (13) 

Lch  Jo 


2.2.  Flow  distribution  model 


The  flow  distribution  model  is  applied  for  both  the  anode  and 
cathode  manifolds.  Mass  and  Momentum  Conservation,  Eqs. 
(14)  and  (18),  are  the  two  main  governing  effects.  While  both 
U-  and  Z-Type  manifolds  are  considered  (see  Fig.  1),  note  that  U- 
Type  manifolds  are  preferred  for  packaging  reasons.  The  anode 


and  cathode  reactants  enter  and  leave  the  inlet  and  outlet  headers 
in  the  same  direction,  for  both  coflow  and  counterflow. 

Mass  conservation 
9  If 

—  (pu)  H - (p  rh  dl  =  0  (14) 

dx  A  JdA 

is  applied  throughout  the  manifolds  and  unit  cells,  where  v  =  y 
along  the  unit  cell  channels  and  v  =  z  along  the  headers.  The 
velocity  u  is  related  to  the  molar  flow  rate  Q  and  the  pressure  P  by 


u  = 


QKT 

PA 


(15) 


where  the  Ideal  Gas  Law  has  been  applied.  In  particular,  Eq. 
(14)  takes  the  form 


q( in)  _  Q(out)  _|_  g(cell) 
q( in)  _j_  Q(ce  11)  _  g(out) 


(16) 

(17) 


at  each  tee  junction  along  the  inlet  and  outlet  header  respectively 
and  for  each  constituent  of  the  gas  stream;  here  g(in)/g(°ut) 
represents  the  flow  entering/exiting  the  tee  junction  along  the 
header,  and  £Tcel6  is  the  unit  cell  flow.  Along  the  unit  cell  chan¬ 
nels,  Eq.  (14)  takes  the  forms  Eqs.  (3)-(6). 

Momentum  conservation 


(p  mu  dl  + 

JdA 


dP 

dx 


/P«2 

2Dh 


(18) 


is  also  applied  throughout  the  manifolds,  where  v  =  y  along  the 
unit  cell  channels  and  v  =  z  along  the  headers.  Here,  §dA  mu 
represents  a  source  or  sink  of  mass  flux  due  to  flow  splitting 
or  recombination  at  header  tee  junctions  and  mass  exchange 
with  the  GDL  along  the  channels.  The  term  fpu2 /(2Z>h)  is  a 
friction  loss  term  which  includes  the  fitted  Darcy  friction  factor 
/  which  can  be  varied  from  cell  channel  to  cell  channel  and 
the  headers.  Assuming  uniform  conditions  along  each  section 
of  channel/header  with  length  Ax  and  bounding  area  A  A,  the 
mass  flux  term  is  given  by 


(p  mu  dl  = 

JdA 


Q 


C  AA 


MQ 2 

C  Ax 


(19) 


where  Ax  is  the  length  of  the  section  being  considered. 


Stack  inlet 


Stack  inlet 


Outlet  Header 

Stack  outlet 


Fig.  1 .  U-Type  manifold  (left)  and  Z-Type  manifold  (right).  The  reactants  are  fed  through  the  stack  inlet  and  distributed  in  parallel  amongst  the  unit  cells.  The  outflow 
from  the  unit  cells  are  combined  in  the  outlet  header  and  exits  through  the  stack  outlet. 
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The  temperature  along  the  inlet  header  and  the  outlet  header 
is  held  constant,  and  taken  to  be  the  same  as  the  coolant  inlet 
and  outlet  values.  A  linear  temperature  profile  is  imposed  along 
the  channels.  Note  that  this  restriction  has  been  removed  in  [13], 
where  the  temperature  distribution  is  computed. 

3.  Computational  strategy 

The  computational  strategy  is  essentially  a  series  of  nested 
iterations.  At  the  outer  level,  Newton  iterations  are  employed  to 
solve  for  the  stoich  distribution  for  both  the  cathode  and  anode 
manifolds.  At  the  inner  level,  the  stoich  distribution  is  fixed  and 
the  header  conditions  and  unit  cell  conditions  are  computed.  The 
molar  flow  rate  and  pressure  along  the  headers  are  computed  us¬ 
ing  Mass  and  Momentum  Conservation,  Eqs.  (14)  and  (18).  The 
unit  cell  conditions,  including  current  density,  voltage,  channel 
pressure,  and  channel  molar  flow  rates,  are  computed  using  the 
unit  cell  model  from  [2] . 

The  unit  cell  channel  outlet  pressure  Pj  is  a  function  of  the 
stoich  distribution  3,  and  this  is  expressed  symbolically  as 

Pj  =  Fj(s ).  (20) 

These  pressures  must  equal  the  pressure  drops  (A P)j  between 
the  jth  and  (j  +  l)th  tee  junctions  along  the  outlet  header  as 
computed  by  Eq.  (18).  The  pressure  drops  (A  P)j  are  computed 
from  a  specified  stoich  distribution  s  by  using  Eqs.  (14)  and  (18) 
to  compute  the  inlet  header  conditions,  using  the  model  from 
[2]  to  compute  the  unit  cell  conditions  (and  in  particular  the 


Table  1 


Comparison  of  the  friction  factors 


/ 

Fitted 

Laminar 

Turbulent 

Exp.  [20] 

Header 

0.207 

0.0004 

0.222 

Channel 

0.117 

0.0543 

0.300 

0.110 

Comparison  of  the  fitted  friction  factor,  the  laminar  friction  factor  computed 
using  Eq.  (33),  the  turbulent  friction  factor  computed  using  Eq.  (34)  with  e  =  0, 
and  the  friction  factor  from  the  experimental  data  from  [20] .  The  fitted  friction 
factors  are  averaged  over  the  length  of  the  channel  and  header  respectively.  The 
laminar  and  turbulent  friction  factors  are  computed  using  stack  inlet  conditions. 

conditions  at  the  unit  cell  channel  outlets),  and  using  Eq.  (18)  to 
compute  (A P)j,  taking  the  pressure  at  the  start  of  each  header 
section  to  be  the  same  as  the  unit  cell  outlet  pressure.  This  gives 
N  —  1  equations  for  the  stoich  distribution  s : 

Gj(s)  =  F j+i(s)  -  Fj(s)  ~  (A P)j(s)  =  0,  (21) 

for  j  =  1 ,  . . . ,  N  —  1 .  In  addition,  the  average  of  the  stoich  dis¬ 
tribution  must  match  the  target  average  stoich  savg: 

N 

Gn(s )  =  y)  sj  -  Nsmg  =  0.  (22) 

7=1 

This  equation  represents  mass  conservation  along  the  header:  the 
sum  of  the  flow  going  into  the  unit  cells  is  equal  to  the  header 
inlet  flow. 

Newton  iterations  are  performed  to  solve  Eqs.  (21)  and  (22) 
for  the  stoich  distribution,  and  takes  the  form 

j(i+D  =  3(0  _  y-i(3(0)G(3(0)>  (23) 
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Fig.  2.  Left:  comparison  of  the  computed  (solid  line)  and  measured  (crosses)  static  pressures  along  the  cathode  inlet  and  outlet  headers  of  a  Mk  7  fuel  cell  stack  for 
current  densities  of  1,  0.88  and  0.77  A  cm“2.  Right:  comparison  of  the  air  flow  distribution  computed  with  the  model  (solid  line)  and  the  empirical  formula  Eq.  (35) 
applied  to  experimental  data.  Other  operating  conditions:  T  =  343  K,  air/H2  stoichiometry  =  2.0/ 1 .5,  dew  point  =  343/347  K  cathode/anode.  For  the  computations, 
the  friction  factors  for  the  header  and  the  flow  channels  were  fitted  to  /h  =  0.207  and  /ch  =  0.1 17  respectively. 
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where  s®  is  the  Zth  iterate  for  the  stoich  distribution,  and  J  is  the 
Jacobian  matrix  of  G  defined  by 


Jij  — 


dGi 

dsj 


(24) 


However,  to  reduce  the  computational  cost  of  computing  /,  it  is 
assumed  that  the  dFj/dsj  are  the  dominant  terms  in  /  and  the 
lower  order  terms  are  dropped  to  obtain  the  approximation 


/_9£l 
/  asi 


dF2 

ds2 

dF2 

ds2 


dF3 

ds3 


\ 


9Fn- 2  dFtf-1 
dsjv-2  dsjv-l 

9Fn-1  dFjy 

dS]\f—i  dSft 

1  •••  1  1  1  / 


(25) 


where  dFj/dsj  is  computed  numerically.  The  resulting  approxi¬ 
mate  Newton  iterations  converge  for  a  wide  range  of  parameters, 
and  are  robust. 

For  the  cathode  manifold,  the  molar  flux  at  the  stack  inlet  for 
oxygen  and  nitrogen  is  taken  as 


Q 


(in) 

o 


I 

4F’ 


/i  (in)  ^  (in) 

*2n  2\  ’ 


(26) 

(27) 


while  for  the  anode  manifold 


Q 


(in) 

h 


I 

2F' 


(28) 


For  both  manifolds,  the  inlet  molar  flux  for  water  is  deter¬ 
mined  by  specifying  the  dew  point  Tdew- 


o(in)  = 

W 


^sat(^dew) 

P  ~  Z>sat(^dew) 


(6oin)  +  en 


(in)' 


(29) 


Q 


(in)  = 
w 


Z)sat(2dew)  ^(in) 


P  ^sat(^dew) 


(30) 


for  the  cathode  and  anode  respectively.  The  unit  cell  channel 
inlet  fluxes  are  then  determined  by  the  stoich  distribution: 


eK  = 


sj 


Ns 


Q? 


avg 


(31) 


where  X  refers  to  oxygen,  nitrogen,  and  water  for  the  cathode, 
and  hydrogen  and  water  for  the  anode. 

For  each  unit  cell,  Newton  iterations  are  employed  to  solve 
for  the  current  density  profile  such  that  the  voltage  is  constant 
in  the  along-the-channel  direction.  Note  that  this  restriction  is 
removed  in  [1,9].  For  numerical  stability  reasons  however,  the 
equation 

d2V 


(32) 


along  with  Eq.  (13)  are  solved  instead.  Employing  a  so¬ 
lution  procedure  analogous  to  that  for  the  stoich  distribu¬ 
tion,  the  current  density  profile  is  thus  solved.  See  [2]  for 
details. 


Polarization  Curve 


Fig.  3.  Polarization  curve  with  varying  cathode  inlet  pressures  of  300,  273.6, 
and  200  kPa.  The  operating  points  corresponding  to  the  three  base  cases  X,  A,  B 
are  indicated  on  the  300  kPa  curve.  Aj  and  Bi  indicate  operating  points  where 
the  cathode  inlet  pressure  for  Case  A  and  B  is  lowered  to  200  and  273.6  kPa 
respectively,  corresponding  to  the  conditions  for  the  runs  shown  in  Figs.  10  and 
11.  All  other  operating  conditions  are  as  listed  in  Table  2. 

The  model  and  solution  procedure  were  implemented  using 
Java  5  and  C.  On  a  AMD64  2.2  GHz  machine  with  1  GB  of 
RAM,  computations  for  a  100  cell  stack  take  under  1  min  to 
complete. 

4.  Experimental  validation 

The  experimental  results  were  obtained  using  a  U- 
manifolded  Ballard  Mk  7  fuel  cell  stack  consisting  of  190 
unit  cells.  Pressure  measurements  were  performed  in  the 
inlet  and  outlet  headers  using  a  custom  built  Pitot-tube  which 
was  inserted  from  the  downstream  end  (relative  to  the  inlet 
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Fig.  4.  Calculated  cell  voltage  vs.  cathode  stoich  for  the  base  cases  A,  B  from 
Fig.  3.  Also  given  are  operating  points  A2  and  B2  corresponding  to  the  reduced 
cathode  stoichs  of  1.13  and  1.5  respectively.  These  operating  conditions  were 
used  for  the  results  shown  in  Figs.  8  and  9.  All  other  conditions  as  in  Table  2. 
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Fig.  5.  Unit  cell  model  runs  corresponding  to  Case  B  which  show  the  effect  of  lowering  the  cathode  stoich  from  1.8  (base  case,  solid  line)  to  1.2  (understoich,  dashed 
line).  Plots  include  (a)  cell  voltage,  (b)  current  density,  (c)  cathode  channel  pressure,  (d)  anode  channel  pressure,  (e)  cathode  channel  velocity,  and  (f)  anode  channel 
velocity.  The  cathode  inlet/outlet  and  anode  outlet/inlet  are  located  at  the  positions  0  m/0.7  m  along  the  channel. 


header)  of  the  stack  for  pressure  measurements  in  the  oxidant 
inlet  and  outlet  headers.  This  arrangement  allowed  for  the 
measurement  of  the  static  and  dynamic  pressures  in  the  inlet 
header  and  the  static  pressure  in  the  outlet  header.  The  cathode 
and  anode  stoichiometries  were  set  to  2.0  and  1.5  respectively, 
and  pressure  measurements  were  performed  at  current  densities 
of  1,  0.88,  and  0.77  A  cm-2.  Fig.  2  shows  the  experimental 
results  together  with  model  calculations  obtained  by  using  the 
friction  factors  in  the  headers  and  flow  channels  as  the  only 
fitting  parameters.  These  parameters  appear  in  Table  1  and  are 
the  same  for  each  current  density  case.  There  is  good  agreement 
between  experimental  measurements  and  model  predictions 
for  the  pressure  variations  in  the  cathode  header.  Losses  due 
to  friction  at  the  tee  junctions  and  bends  are  neglected  as 
these  losses  are  insignificant  for  Re  >  100  [14].  This  regime 
Re  >100  corresponds  to  a  current  density  of  at  least  0.15  and 
0.33  A  cm-2  for  the  oxidant  (1.8  stoichiometry)  and  fuel  (1.2 
stoichiometry)  streams  respectively.  These  current  densities  are 
significantly  smaller  than  those  at  which  flow  distribution  has 
an  impact  on  cell  performance  (>  0.5  A  cm  ,  see  Section  5). 

Table  1  compares  the  fitted  friction  factors  for  the  headers  and 
channels  versus  that  obtained  using  the  laminar  flow  equation 


(33) 


versus  the  Colebrook  equation  for  turbulent  flow 


1  (e/Dh  ,  2.51  \ 

77  V3.7+^v?j’ 


(34) 


where  e/D^is  the  relative  roughness  of  the  pipe.  The  friction  fac¬ 
tors  were  obtained  at  the  inlet  of  the  header  and  channel  lengths 


using  the  actual  channel  and  header  geometry.  The  fitted  fric¬ 
tion  factors  are  between  the  laminar  and  turbulent  values  ob¬ 
tained  from  Eqs.  (33)  to  (34),  indicating  that  the  fitted  friction 
factors  are  reasonable.  It  is  seen  that  the  flow  is  in  the  mixed 
zone  in  the  channels,  and  turbulent  in  the  headers.  In  reality,  the 
friction  factor  varies  along  the  channel  due  to  varying  condi¬ 
tions  along  the  channel  and  the  appearance  of  two-phase  flow; 

Table  2 


Base  case  parameters  and  operating  conditions 


Symbol 

Description 

Value 

Ach 

Cross-sectional  area  of  channel  (m2) 

3.00  x  10-7 

Ah 

Cross-sectional  area  of  header  (m2) 

3.24  x  10“4 

/ch 

Friction  factor  of  channel  (— ) 

0.1 

A 

Friction  factor  of  header  (— ) 

0.2 

i 

Average  current  density  (Case  X/A/B)  (A  cm-2) 

0.10/0.92/1.33 

Tch 

Length  of  channel  (m) 

0.7 

Lh 

Length  of  header  (m) 

0.4 

N 

Number  of  unit  cells  (— ) 

100 

Nch 

Number  of  channels  per  unit  cell  (— ) 

36 

Pch 

Cross-sectional  perimeter  of  channel  (m) 

2.19  x  10“3 

Ph 

Cross-sectional  perimeter  of  header  (m) 

7.20  x  10“2 

p(  in) 

Inlet  pressure  (anode/cathode)  (kPa) 

320/300 

Average  stoichiometry  (anode/cathode)  (— ) 

1.2/1. 8 

T 

Temperature  (inlet/outlet)  (K) 

343/353 

Tdew 

Dew  point  at  inlet  (anode/cathode)  (K) 

347/343 

— 

Flow  direction  of  anode  vs.  cathode  gases 

Counterflow 

— 

Manifold  type 

U-Type 

Base  case  parameters  and  operating  conditions  for  the  parameter  sensitivity 
analysis.  Air  and  H2  are  used  as  the  oxidant  and  fuel  respectively,  and  their  tem¬ 
perature  profiles  are  taken  to  be  the  same  as  the  coolant.  The  coolant  temperature 
is  held  constant  along  the  headers,  and  varies  linearly  along  the  channel  direc¬ 
tion.  The  coolant  flow  is  coflow  with  respect  to  the  oxidant,  and  counterflow 
with  respect  to  the  fuel.  The  cathode  inlet  is  located  at  0  m. 
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Fig.  6.  (a)  Calculated  flow  distribution,  (b)  average  voltage  distributions,  (c)  current  density  distribution,  and  (d)  cathode  header  pressures  as  functions  of  cell  number 
for  a  U-Type  manifolded  stack  for  Case  A  (see  Figs.  3  and  4).  In/outlets  at  cell  #1. 


however,  the  use  of  the  averaged  friction  factor  is  a  good  first 
approximation. 

Currently,  a  direct  method  to  measure  the  flow  distribution 
does  not  appear  to  have  been  demonstrated.  As  a  result,  only 
indirect  measurements  are  available  requiring  the  use  of  approx¬ 
imate  equations  such  as 


^avg 


to  derive  the  actual  flow  distribution.  Eq.(35)  is  derived  from 


DhAP/  Ax 
pu2 /2 


(36) 


by  assuming  f  D h,  Ax,  and  p  are  the  same  for  all  channels  and 
cells.  Direct  pressure  measurements  in  the  header  are  notably 
difficult  to  perform  because  of  interferences  from  two-phase 
flow  (especially  at  the  outlet  header),  the  requirement  for  rela¬ 
tively  long  stacks  (may  require  other  stack  design  changes  with 
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respect  to  cell  alignment  and  compression)  and  sensitive  pres¬ 
sure  sensors.  These  last  requirements  are  necessary  since  headers 
are  designed  such  that  their  pressure  drop  is  significantly  smaller 
than  the  cell  flow  field  pressure  drop  to  ensure  a  uniform  flow 
distribution.  Fig.  2  shows  good  agreement  between  the  flow  dis¬ 
tribution  calculated  with  the  empirical  formula  Eq.  (35)  and  the 
model  solution. 

Independent  measurements  of  oxidant  pressure  drop  versus 
average  gas  velocity  were  also  obtained  for  the  Ballard  Mk  7 


stack  [20].  These  data  were  correlated  to  a  second  order  poly¬ 
nomial  of  the  form 

AP  =  1725.3m  +  38.9m2;  (37) 

The  corresponding  friction  factor  is  estimated  by  substituting 
Eq.  (37)  into  Eq.  (36),  yielding 

/  =  4.4488  x  10-4(1725.3m-1  +  38.9).  (38) 
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For  a  cathode  stoichiometry  of  2,  the  average  gas  velocity 
is  8.34  ms-1,  and  so  the  friction  factor  is  0.110.  This  estimate 
further  indicates  that  the  fitted  friction  factor  for  the  channel  is 
reasonable  (see  Table  1). 

5.  Sample  runs 

Table  2  lists  the  base  case  parameters  and  operating  condi¬ 
tions  which  shall  be  considered  for  all  sample  runs.  Three  base 
case  conditions  shall  be  considered:  Case  X  in  the  activation 
region  with  i  =  0.10  A  cm-2,  Case  A  in  the  ohmic  region  with 
i  =  0.92  A  cm-2  and  Case  B  near  the  mass  transport  limited  re¬ 


gion  with  i  =  1 .33  A  cm-2 .  The  polarization  curve  generated  by 
the  unit  cell  model  is  shown  in  Fig.  3,  and  points  correspond¬ 
ing  to  these  base  cases  are  indicated.  In  addition,  the  polarization 
curve  corresponding  to  a  reduced  cathode  inlet  pressure  of  273.6 
and  200  kPa  are  shown.  In  Fig.  4,  voltage  versus  cathode  stoi¬ 
chiometry  graphs  for  a  single  unit  cell  are  shown.  These  graphs 
exhibit  a  typical  threshold  behaviour  in  which  the  minimum  sto¬ 
ich  needed  by  a  typical  unit  cell  increases  with  current  density. 

In  Fig.  5,  the  unit  cell  model  was  used  to  investigate  the  ef¬ 
fect  of  lowering  the  cathode  stoich  from  1.8  to  1.2  for  Case 
B.  As  the  cathode  flow  decreases,  the  performance  of  the  unit 
cell  suffers  due  to  the  increasing  scarcity  of  oxygen  (Fig.  5a), 
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and  the  current  density  distribution  shifts  towards  the  cath¬ 
ode  inlet  at  0  m  (Fig.  5b).  The  decreased  cathode  channel  in¬ 
let  flow  leads  to  a  lower  pressure  drop  and  velocity  as  ex¬ 
pected  (Fig.  5c  and  e).  However,  because  the  hydrogen  is  flow¬ 
ing  counter  to  the  oxygen,  more  hydrogen  is  consumed  near 
the  anode  outlet,  leading  to  more  flow  through  the  anode  chan¬ 
nel  and  therefore  increased  anode  pressure  drop  (Fig.  5d  and 
f).  There  is  little  along-the-channel  variation  in  the  cathode 
channel  velocity  due  to  the  large  presence  of  non-reactive  ni¬ 
trogen  (Fig.  5e),  but  the  anode  channel  velocity  decreases  sig¬ 


nificantly  from  anode  inlet  to  outlet  as  most  of  the  hydrogen  is 
consumed. 

In  discussing  the  following  sample  runs,  it  is  helpful  to  sum¬ 
marize  the  following  principles  which  govern  steady  state  flow 
distribution: 

1.  The  pressure  drop  of  the  gases  flowing  from  stack  inlet  to 
outlet  must  be  the  same  along  all  paths  taken. 

2.  The  pressure  drop  of  the  gases  increases  with  increasing  flow 
rate,  friction  factor,  and  distance  travelled. 
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3.  Decreasing  cathode  flow  in  a  unit  cell  leads  to  a  decreased 
cathode  pressure  drop  across  that  cell,  but  an  increased 
anode  pressure  drop,  under  counterflow  operation  (see 
Fig.  5c  and  d). 

Base  case  results  for  Cases  A  and  B  are  shown  in  Figs.  6 
and  7.  Because  the  manifolding  is  U-Type,  the  stack  inlet  and 
outlet  are  located  close  to  the  first  unit  cell.  The  cathode  flow 
distribution  exhibits  a  greater  variation  than  the  anode,  due  to 


the  fluid  properties  and  reactant  flow  rates.  Geometry,  such  as 
channel  cross-sectional  area  or  length,  plays  no  role  in  this  study 
as  the  anode  and  cathode  geometry  are  the  same.  Both  the  anode 
and  cathode  flow  distributions  tend  to  favor  the  unit  cells  which 
are  closer  to  the  stack  inlet,  because  of  the  shorter  travel  distance. 
However,  the  Bernoulli  effect  shifts  the  anode  flow  distribution 
to  the  right,  since  the  increasing  static  pressure  along  the  inlet 
header  induces  increasing  pressure  drops  across  the  unit  cells, 
thus  drawing  more  flow.  As  a  result,  these  two  competing  effects 
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result  in  a  minimum  in  the  anode  flow  distribution  around  cell 
#20 — the  ‘shorter  distance’  effect  appears  to  have  more  impact 
on  the  anode  flow  near  stack  inlet,  while  the  Bernoulli  effect 
dominates  farther  from  the  stack  inlet. 

While  the  cell  voltage  is  virtually  flat  for  Cases  X  (not  shown) 
and  A,  there  is  somewhat  more  variation  for  Case  B  as  that  op¬ 
erating  point  is  closer  to  the  mass  transport  knee  in  the  polar¬ 
ization  curve.  The  variation  of  current  density  along  the  channel 
direction  is  smaller  for  cells  closer  to  the  stack  inlet  due  to  the 
increased  flow.  The  static  pressure  along  the  cathode  inlet  header 
increases  slightly  because  of,  again,  the  Bernoulli  effect. 
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In  Figs.  8  and  9,  the  cathode  stoichs  for  Cases  A  and  B  are 
lowered  to  1.13  and  1.5  respectively  (Cases  A2  and  B2).  This  is 
equivalent  to  moving  from  the  ohmic  regime  to  the  mass  trans¬ 
port  limited  regime.  In  both  cases,  these  operating  conditions  are 
closer  to  the  knee  in  Fig.  4.  Although  the  flow  distributions  do  not 
change  appreciably,  the  average  cell  voltage  distribution  in  the 
stack  changes  drastically:  the  variation  in  average  cell  voltage 
increases  from  3  to  70  mV  and  9  to  60  mV  respectively.  The  cell 
voltage  variation  is  further  exacerbated  if  manufacturing  varia¬ 
tions  are  taken  into  account,  as  will  be  discussed  later.  For  the 
base  cases,  the  cell  voltage  decreases  monotonically  with  cell 


Anode  Flow 


x  <]  o4Current  Density  (Cell  50) 


Fig.  17.  Effect  of  variation  of  the  cathode  channel  friction  factor  for  the  inlet  and  outlet  headers  on  operating  point  B.  Text:  see  Fig.  14. 
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Cathode  Flow  Anode  Flow 
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Fig.  18.  Effect  of  a  random  variation  (20%  S.D.)  of  the  cathode  channel  friction  factor  on:  (a)  cathode  flow  distribution,  (b)  anode  flow  distribution,  (c)  average  cell 
voltage,  and  (d)  average  current  density  distribution  for  operating  point  X  (see  Fig.  3).  In  (d),  the  cathode  inlet/outlet  is  at  the  positions  0  m/0.7  m  along  the  channel. 


number,  whereas  for  Cases  A2  and  B2,  the  cell  voltage  shows  a 
shallow  minimum  around  cell  #80  (see  Figs.  6b,  7b,  8b  9b). 

In  Figs.  10  and  11,  the  cathode  inlet  pressure  is  lowered  to 
200  and  273.6  kPa  for  Cases  Ai  and  Bi  respectively  (see  Fig.  3). 
There  is  a  greater  pressure  drop  along  the  channels  since  the  gas 
velocity  increases  in  the  manifold  for  these  lower  inlet  pressures. 
As  in  the  previous  case  with  reduced  stoichiometry,  the  variation 
in  cell  voltage  increases  drastically  to  140  and  140  mV  for  Cases 
Ai  and  Bi  respectively.  It  is  interesting  to  note  that,  for  Case  B 
(the  higher  current  density  case),  the  cell  voltage  distribution 


appears  to  be  much  more  sensitive  to  pressure  than  for  Case 
A.  This  is  indicated  by  the  fact  that  a  25  kPa  reduction  in  Case 
B  has  a  comparable  effect  as  100  kPa  in  Case  A.  Therefore, 
control  of  operating  conditions  is  critical  to  achieving  optimal 
performance,  especially  at  high  current  densities  located  in  the 
mass  transport  limited  regime. 

In  Fig.  12,  the  manifold  type  for  Case  B  is  changed  from  U- 
Type  to  Z-Type.  The  stack  inlet  is  located  close  to  the  first  cell, 
and  the  stack  outlet  is  located  close  to  the  last  cell.  Anode  and 
cathode  inlet  and  outlet  gases  all  flow  in  the  same  direction.  The 
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Fig.  19.  Effect  of  a  random  variation  (20%  S.D.)  of  the  cathode  channel  friction  factor  on  operating  point  A.  Text:  see  Fig.  18. 
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Fig.  20.  Effect  of  a  random  variation  (20%  S.D.)  of  the  cathode  channel  friction  factor  on  operating  point  B.  Text:  see  Fig.  18. 


anode  and  cathode  flow  distributions  favor  the  unit  cells  which 
are  closer  to  the  stack  outlet.  As  all  gases  travel  the  same  distance 
from  stack  inlet  to  outlet  along  all  paths,  the  ‘shorter  distance’ 
effect  does  not  apply.  Instead,  because  of  the  large  pressure  drop 
along  the  outlet  header  due  to  the  accumulation  of  unused  re¬ 
actant  and  product  gases  (see  Figs.  12d  and  13d,  for  instance), 
the  pressure  drop  across  the  last  unit  cell  must  be  quite  large  in 
order  that  the  steady  state  pressure  distributions  be  consistent. 
Consequently,  more  flow  must  be  directed  to  unit  cells  closer 
to  the  stack  outlet.  There  is  a  greater  spread  in  the  flow  distri¬ 
bution  compared  to  the  U-Type  manifold,  and  this  is  reflected 
in  the  somewhat  greater  spread  in  cell  voltage.  In  Fig.  13,  the 

Cathode  Flow 
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cathode  inlet  pressure  is  lowered  to  273.6  kPa.  Compared  with 
the  U-Type  manifold,  the  performance  of  the  Z-Type  manifold 
is  seen  to  be  worse  -  there  is  a  140  mV  spread  in  cell  voltage  for 
the  U-Type  manifold  versus  330  mV  for  the  Z-Type  manifold. 
Therefore,  U-Type  manifolds  are  preferable  on  a  performance 
basis,  in  addition  to  the  packaging  basis.  It  is  interesting  to  note 
the  change  in  current  density  distribution  along  the  length  of 
the  stack.  Whereas  for  the  U-Type  manifold,  the  current  den¬ 
sity  at  the  reactant  inlets/outlets  shows  increasing/decreasing 
behaviour  with  the  cell  number,  this  trend  is  reversed  for  the 
Z-Type  manifold.  But,  as  for  the  U-Type  manifold,  the  high¬ 
est  cell  voltage  is  still  associated  with  the  highest  oxidant  flows 
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Fig.  21.  Effect  of  a  random  variation  (20%  S.D.)  of  the  cathode  channel  friction  factor  on  operating  point  Bi.  Text:  see  Fig.  18. 
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and  with  the  lowest  current  density  variation  inlet  — >►  middle  -> 
outlet  within  the  cells  (see  Figs,  lla-c  and  13a-c). 

In  Figs.  14  and  15,  the  cathode  channel  friction  factor  for  the 
center  cell  is  varied  between  0.05,  0.10  (base  case),  and  0.15. 
The  anode  channel  friction  factors  were  fixed  at  the  base  case 
value.  The  anode  flow  is  affected  by  the  changing  cathode  fric¬ 
tion  factors  in  the  same  way  as  on  the  cathode,  albeit  to  a  lesser 
extent.  This  is  due  to  the  coupling  of  the  cathode  and  anode 
gas  flows  through  reactant  consumption  and  electrical  current. 
This  is  an  important  finding  that  shows  that  flow  distribution  in 
a  given  electrode  compartment  cannot  be  considered  on  its  own, 
but  in  combination  with  the  other  compartment  because  of  cou¬ 
pling  between  the  reactant  compartments.  This  is  an  effect  that 
has  apparently  not  been  previously  noted  or  modelled.  While 
the  cell  voltage  for  Case  A  is  not  affected  much  by  the  chang¬ 
ing  friction  factor,  Case  B  is  to  a  much  larger  extent.  Higher 
channel  resistance  also  leads  to  a  shift  of  the  current  density 
profile  towards  the  channel  inlet.  The  differing  current  density 
profile  between  the  anomalous  cell  and  its  neighbors  leads  to  in¬ 
plane  currents  in  the  bipolar  plates  which  the  model  presented 
here  does  not  capture.  However,  this  important  effect  has  been 
incorporated  in  other  models  [1,9]. 

In  Fig.  16,  the  cathode  channel  friction  factor  for  all  cells 
is  varied  for  Case  B.  While  higher  channel  resistances  lead  to  a 
more  uniform  flow  distribution,  this  would  incur  more  losses  due 
to  higher  reactant  pressure  losses.  In  Fig.  17,  the  cathode  header 
friction  factor  is  varied  for  Case  B.  Lower  header  resistances 
lead  to  a  more  uniform  flow  distribution.  These  runs  illustrate 
the  use  of  the  model  to  quantitatively  balance  requirements  for 
a  series  of  different  scenarios,  such  as  attaining  a  uniform  flow 
distribution  versus  reducing  compressive  power. 

Normally  distributed  random  variations  in  the  friction  factor 
in  the  cathode  channels  are  introduced  with  20%  S.D.  for  all  base 
cases,  and  the  results  are  shown  in  Figs.  18-20.  These  variations 
might  be  the  result  of  manufacturing  process  control  limitations. 
The  cathode  flow  distribution  is  virtually  the  same  in  all  three 
cases,  but  the  anode  flow  distribution  is  more  spread  out  for  Case 
B  than  for  Cases  X  and  A.  The  cathode  flow  distribution  induces 
a  cell  voltage  distribution  with  a  spread  of  around  0.008  V  for 
Case  X,  0.015  V  for  Case  A,  and  0.4  V  (!)  for  Case  B.  These 
results  indicate  that,  while  variations  in  cathode  flow  distribu¬ 
tion  have  minor  effect  in  the  kinetic  and  ohmic  regions  of  the 
polarization  curve,  these  variations  can  have  a  huge  effect  in  the 
mass  transport  limited  region.  The  existence  of  a  low  cell  volt¬ 
age  in  a  single  cell  is  sufficient  to  create  a  fault  condition  which 
could  lead  to  a  stack  shutdown.  In  the  case  of  the  cathode,  harm¬ 
ful  effects  are  less  than  on  the  anode  [19],  but  shutdown  cannot 
be  prevented  unless  other  precautions  are  taken  to  identify  the 
origins  of  the  cell  failure  (detection  of  hydrogen  in  the  cathode 
exhaust,  for  example). 

Fig.  21  illustrates  the  combined  effect  of  normally  distributed 
random  manufacturing  variations  in  cathode  channel  friction 
factor  with  a  reduction  in  cathode  compartment  pressure  (a 
more  realistic  and  practical  case).  There  is  a  significant  in¬ 
crease  in  cell  performance  variation,  from  0.4  V  (Fig.  20c)  to 
2  V  (Fig.  21c),  even  though  the  cathode  and  anode  flow  dis¬ 
tributions  have  hardly  changed  (Figs.  20a  and  b  and  21a  and 


b).  A  significant  increase  in  current  density  distribution  is  also 

_ 

observed,  with  a  spread  increasing  from  0.8  to  1.2  A  cm 
(Figs.  20d  and  2 Id).  The  negative  cell  voltages  observed  are 
an  artifact  of  the  present  model  because  the  model  does  not 
account  for  additional  electrochemical  reactions  which  would 
limit  the  cell  voltage  to  approximately  0  V  by  the  appearance 
of  hydrogen  evolution  [19].  Fig.  21  shows  that  control  of  op¬ 
erating  conditions  and  manufacturing  processes  is  important  to 
ensure  uniform  performance  in  a  stack,  especially  if  design  at 
high  current  densities  is  pursued  and  other  design  (membrane 
thickness,  catalyst  loading,  etc.)  and  operating  conditions  (re¬ 
actant  flow  rate,  cell  temperature,  etc.)  variations  are  taken  into 
account. 

6.  Conclusion 

A  PEMFC  stack  model  which  couples  flow  distribution  with 
unit  cell  performance  has  been  presented.  Mass  and  momentum 
conservation  were  applied  to  solve  for  the  flow  and  pressure 
distribution  of  the  fuel  and  oxidant  manifolds,  while  an  earlier 
unit  cell  model  developed  by  our  group  [2]  was  used  to  de¬ 
termine  cell  performance.  Individual  parameters  for  each  unit 
cell  can  be  varied,  and  so  the  model  can  be  used  to  investi¬ 
gate  the  effects  of  different  designs,  operating  conditions,  and 
manufacturing  requirements  to  minimize  the  impact  on  stack 
operability.  For  instance,  the  model  can  be  used  to  identify 
regions  of  high  current  density  that  could  induce  stack  fail¬ 
ure  or  limit  life  by  exceeding  the  membrane  glass  transition 
temperature. 

A  parameter  sensitivity  analysis  was  also  performed  in  which 
the  header  and  channel  friction  factors  were  varied,  and  the  effect 
on  performance  was  determined.  In  general,  higher  channel  fric¬ 
tion  factors  lead  to  a  more  uniform  flow  distribution,  consistent 
with  the  experimental  formula  Eq.  (35),  but  this  in  turn  requires 
more  reactant  compression.  The  coupling  of  the  anode  flow  dis¬ 
tribution  to  the  cathode  via  the  ME  A  was  demonstrated;  this 
effect  has  apparently  not  been  previously  reported  or  observed. 
The  performance  of  the  Z-Type  manifold  was  also  seen  to  be 
worse  than  the  U-Type  manifold.  At  higher  current  densities, 
performance  is  more  sensitive  to  operating  conditions,  such  as 
cathode  stoich  and  inlet  pressure.  This  was  especially  true  when 
random  variations  in  the  cathode  friction  factor  were  introduced 
to  simulate  manufacturing  variability,  where  the  performance  of 
a  few  unit  cells  was  seen  to  suffer  tremendously.  This  highlights 
a  possible  mechanism  wherein  failure  in  a  few  cells  can  lead  to 
failure  of  the  entire  stack. 

Significant  additional  work  is  required  for  model  implemen¬ 
tation.  Most  importantly,  validation  data  linking  manufactur¬ 
ing  parameter  distributions  to  associated  performance  data  is 
needed  to  establish  useful  relationships  for  stack  design  and 
manufacturing.  This  work  will  largely  benefit  from  the  devel¬ 
opment  of  a  direct  method  to  measure  stack  flow  distribution. 
Additional  work  is  planned  to  establish  valid  two-phase  flow 
correlations  to  capture  an  aspect  that  was  disregarded.  This  can 
be  achieved  by  correlating  pressure  drop  and  liquid  water  con¬ 
tent  (image  analysis)  measurements  obtained  by  visualization 
techniques  [6,16,18,21]  for  different  operating  conditions  and 
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cell  designs.  Hydrogen  fuel  re-circulation  is  another  considera¬ 
tion  that  requires  model  modifications.  Finally,  there  is  evidence 
that  although  the  study  of  individual  cell  interactions  is  valuable, 
a  stack  model  including  all  interactions  (including  electrical  and 
thermal  coupling)  is  required  as  experimentally  observed  by  Burt 
et  al.  [5]  and  in  the  present  work. 
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